Five different contrasts were performed in order to understand the effects of cerulenin in C. reinhardtii cells. Those were:
Effect of cerulenin in wild-type (0h vs 4h)
Gene obtention
resWT = results(dds, contrast=c("condition","4h","0h"), cooksCutoff = FALSE,
independentFiltering = FALSE)
resWT
## log2 fold change (MLE): condition 4h vs 0h
## Wald test p-value: condition 4h vs 0h
## DataFrame with 13977 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue
## <numeric> <numeric> <numeric> <numeric> <numeric>
## Cre14.g619950.v5.5 307.8280 0.0281565 0.717068 0.0392662 0.968678
## Cre01.g048900.v5.5 76.4131 0.3799222 1.828657 0.2077602 0.835416
## Cre12.g543400.v5.5 428.9410 -0.2585408 0.320513 -0.8066479 0.419869
## Cre16.g685837.v5.5 255.7288 0.2418724 0.453962 0.5328032 0.594170
## Cre06.g278159.v5.5 9770.0557 0.1504420 0.124149 1.2117907 0.225592
## ... ... ... ... ... ...
## Cre06.g278104.v5.5 889.54623 -0.4304988 0.309745 -1.3898503 0.1645743
## Cre17.g717150.v5.5 146.81853 0.0579478 1.170897 0.0494901 0.9605287
## Cre10.g443600.v5.5 170.10265 -1.9623294 1.081889 -1.8137987 0.0697087
## Cre14.g621400.v5.5 5.13639 -0.1021788 4.276178 -0.0238949 0.9809364
## Cre10.g441750.v5.5 1599.45474 0.1400955 0.181245 0.7729615 0.4395452
## padj
## <numeric>
## Cre14.g619950.v5.5 0.997405
## Cre01.g048900.v5.5 0.990087
## Cre12.g543400.v5.5 0.862068
## Cre16.g685837.v5.5 0.942215
## Cre06.g278159.v5.5 0.652547
## ... ...
## Cre06.g278104.v5.5 0.550696
## Cre17.g717150.v5.5 0.996460
## Cre10.g443600.v5.5 0.324773
## Cre14.g621400.v5.5 0.997405
## Cre10.g441750.v5.5 0.870623
summary(resWT)
##
## out of 13977 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up) : 909, 6.5%
## LFC < 0 (down) : 980, 7%
## outliers [1] : 0, 0%
## low counts [2] : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
log.fold.change_WT <- resWT$log2FoldChange
q.value_WT <- resWT$padj
names(log.fold.change_WT) <- gene.ids
names(q.value_WT) <- gene.ids
base.meanWT <- resWT$baseMean
names(base.meanWT) <- gene.ids
# Take a look at upregulated and downregulated genes:
top_20_overexpressed_genes_WT <- resWT[order(resWT$log2FoldChange, decreasing = T),]
kable(top_20_overexpressed_genes_WT[1:20,-(3:4)])
|
|
baseMean
|
log2FoldChange
|
pvalue
|
padj
|
|
Cre10.g465763.v5.5
|
17.445204
|
20.612283
|
0.0000012
|
0.0000297
|
|
Cre07.g331350.v5.5
|
19.178658
|
20.226262
|
0.0000004
|
0.0000097
|
|
Cre10.g466200.v5.5
|
15.235117
|
18.706363
|
0.0000109
|
0.0002224
|
|
Cre02.g143107.v5.5
|
7.114070
|
17.784866
|
0.0000295
|
0.0005527
|
|
Cre06.g250650.v5.5
|
13.336058
|
15.503286
|
0.0002575
|
0.0037533
|
|
Cre02.g102300.v5.5
|
33.165537
|
8.911089
|
0.0019320
|
0.0205662
|
|
Cre17.g735876.v5.5
|
22.382807
|
8.046166
|
0.0126038
|
0.0945220
|
|
Cre02.g116900.v5.5
|
12.062879
|
7.629730
|
0.0605932
|
0.2956062
|
|
Cre08.g364850.v5.5
|
8.610772
|
7.489975
|
0.0320338
|
0.1896386
|
|
Cre16.g680454.v5.5
|
350.851311
|
7.410076
|
0.0001349
|
0.0021043
|
|
Cre02.g091226.v5.5
|
16.082502
|
7.296935
|
0.0859688
|
0.3740927
|
|
Cre12.g511952.v5.5
|
11.269179
|
7.184714
|
0.0909203
|
0.3867668
|
|
Cre02.g095151.v5.5
|
13189.468430
|
7.063199
|
0.0000000
|
0.0000000
|
|
Cre07.g341750.v5.5
|
36.286173
|
6.880834
|
0.0077292
|
0.0644578
|
|
Cre08.g369200.v5.5
|
5.101397
|
6.578353
|
0.1158762
|
0.4496394
|
|
Cre12.g540351.v5.5
|
46.825918
|
6.539815
|
0.0033205
|
0.0323643
|
|
Cre03.g167690.v5.5
|
3.495144
|
6.501141
|
0.0738629
|
0.3369917
|
|
Cre03.g144444.v5.5
|
8.399966
|
6.444530
|
0.1295947
|
0.4817355
|
|
Cre04.g229422.v5.5
|
3.341127
|
6.298153
|
0.0952044
|
0.3978966
|
|
Cre12.g554929.v5.5
|
123.196764
|
6.295387
|
0.0000011
|
0.0000260
|
top_20_overrepressed_genes_WT <- resWT[order(resWT$log2FoldChange, decreasing = F),]
kable(top_20_overrepressed_genes_WT[1:20,-(3:4)])
|
|
baseMean
|
log2FoldChange
|
pvalue
|
padj
|
|
Cre01.g043250.v5.5
|
102.22785
|
-22.99906
|
0e+00
|
0.0e+00
|
|
Cre16.g658300.v5.5
|
46.21710
|
-22.89138
|
0e+00
|
0.0e+00
|
|
Cre02.g115400.v5.5
|
20.74418
|
-22.35411
|
0e+00
|
9.0e-07
|
|
Cre12.g500250.v5.5
|
28.14109
|
-22.03587
|
0e+00
|
0.0e+00
|
|
Cre03.g179400.v5.5
|
39.29917
|
-22.00429
|
0e+00
|
0.0e+00
|
|
Cre06.g292249.v5.5
|
18.96145
|
-21.81929
|
1e-07
|
2.8e-06
|
|
Cre06.g283000.v5.5
|
20.54063
|
-21.80617
|
0e+00
|
1.0e-07
|
|
Cre17.g735750.v5.5
|
27.85918
|
-21.60836
|
0e+00
|
0.0e+00
|
|
Cre01.g034380.v5.5
|
38.57187
|
-21.42391
|
0e+00
|
0.0e+00
|
|
Cre13.g588250.v5.5
|
49.57133
|
-21.38082
|
0e+00
|
0.0e+00
|
|
Cre16.g684715.v5.5
|
14.85207
|
-21.37278
|
4e-07
|
9.9e-06
|
|
Cre03.g143827.v5.5
|
66.50156
|
-21.31874
|
0e+00
|
0.0e+00
|
|
Cre14.g629840.v5.5
|
37.01919
|
-21.31223
|
0e+00
|
0.0e+00
|
|
Cre14.g632350.v5.5
|
93.48560
|
-21.17466
|
0e+00
|
0.0e+00
|
|
Cre12.g554100.v5.5
|
20.57133
|
-21.16066
|
0e+00
|
2.0e-07
|
|
Cre09.g402050.v5.5
|
54.88133
|
-21.12392
|
0e+00
|
0.0e+00
|
|
Cre02.g083050.v5.5
|
18.89988
|
-21.07576
|
1e-07
|
4.0e-06
|
|
Cre09.g397050.v5.5
|
47.73692
|
-20.98995
|
0e+00
|
0.0e+00
|
|
Cre16.g658400.v5.5
|
59.75973
|
-20.89439
|
0e+00
|
0.0e+00
|
|
Cre04.g226550.v5.5
|
12.19448
|
-20.86036
|
9e-07
|
2.3e-05
|
# Differentialy activated and repressed genes with a log fold change and p-adj threshold of 1.5 and 0.05
activated.genes.deseq2_WT <- gene.ids[log.fold.change_WT > 1 & q.value_WT < 0.05]
activated.genes.deseq2_WT <- activated.genes.deseq2_WT[!is.na(activated.genes.deseq2_WT)]
activated.genes.deseq2_WT_good <- substr(activated.genes.deseq2_WT,1,nchar(activated.genes.deseq2_WT)-5)
write.table(activated.genes.deseq2_WT_good,file="activated.genes.deseq2_WT.txt",sep="\t",quote=F,row.names = F)
repressed.genes.deseq2_WT <- gene.ids[log.fold.change_WT < - 1 & q.value_WT < 0.05]
repressed.genes.deseq2_WT <- repressed.genes.deseq2_WT[!is.na(repressed.genes.deseq2_WT)]
repressed.genes.deseq2_WT_good <- substr(repressed.genes.deseq2_WT,1,nchar(repressed.genes.deseq2_WT)-5)
write.table(repressed.genes.deseq2_WT_good,file="repressed.genes.deseq2_WT.txt",sep="\t",quote=F,row.names = F)
length(activated.genes.deseq2_WT_good)
## [1] 458
length(repressed.genes.deseq2_WT_good)
## [1] 350
MA plot
plotMA(resWT,alpha=0.05, ylim=c(-30,30), colSig="grey")
abline(h=c(-1,1), col="black", lwd=1)
points(x = base.meanWT[activated.genes.deseq2_WT],
y = log.fold.change_WT[activated.genes.deseq2_WT],col="red",cex=0.2,pch=19)
points(x = base.meanWT[repressed.genes.deseq2_WT],
y = log.fold.change_WT[repressed.genes.deseq2_WT],col="blue",cex=0.2,pch=19)
text(x =100000, y = 15, label = length(activated.genes.deseq2_WT_good),
col = "red", # Color del texto
font = 2, # Estilo
cex = 1.5) # Tamaño
text(x =100000, y = -15, label = length(repressed.genes.deseq2_WT_good),
col = "blue", # Color del texto
font = 2, # Estilo
cex = 1.5) # Tamaño
legend("bottomright", legend = c("Overexpressed", "Repressed"),
lwd = 3, col = c("red", "blue"))

Effect of cerulenin in mutant (0h vs 4h)
Gene obtention
res_atg8 <- results(dds, list(c("genotypeatg8.condition4h","condition_4h_vs_0h")),
cooksCutoff = FALSE, independentFiltering = FALSE)
res_atg8
## log2 fold change (MLE): genotypeatg8.condition4h+condition_4h_vs_0h effect
## Wald test p-value: genotypeatg8.condition4h+condition_4h_vs_0h effect
## DataFrame with 13977 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue
## <numeric> <numeric> <numeric> <numeric> <numeric>
## Cre14.g619950.v5.5 307.8280 -0.0177850 0.721654 -0.0246447 9.80338e-01
## Cre01.g048900.v5.5 76.4131 -0.3532859 1.831528 -0.1928914 8.47044e-01
## Cre12.g543400.v5.5 428.9410 -0.0527444 0.321006 -0.1643096 8.69487e-01
## Cre16.g685837.v5.5 255.7288 0.3030096 0.456488 0.6637846 5.06828e-01
## Cre06.g278159.v5.5 9770.0557 0.5915680 0.124397 4.7554709 1.97984e-06
## ... ... ... ... ... ...
## Cre06.g278104.v5.5 889.54623 0.402382 0.311428 1.292054 0.196338
## Cre17.g717150.v5.5 146.81853 -1.580214 1.173164 -1.346968 0.177991
## Cre10.g443600.v5.5 170.10265 -0.202010 1.079760 -0.187088 0.851592
## Cre14.g621400.v5.5 5.13639 -2.283619 3.977653 -0.574112 0.565892
## Cre10.g441750.v5.5 1599.45474 0.148161 0.181933 0.814375 0.415430
## padj
## <numeric>
## Cre14.g619950.v5.5 9.98339e-01
## Cre01.g048900.v5.5 9.70421e-01
## Cre12.g543400.v5.5 9.75746e-01
## Cre16.g685837.v5.5 8.13498e-01
## Cre06.g278159.v5.5 3.98446e-05
## ... ...
## Cre06.g278104.v5.5 0.517192
## Cre17.g717150.v5.5 0.490686
## Cre10.g443600.v5.5 0.972125
## Cre14.g621400.v5.5 0.847746
## Cre10.g441750.v5.5 0.747101
log.fold.change_atg8 <- res_atg8$log2FoldChange
q.value_atg8 <- res_atg8$padj
names(log.fold.change_atg8) <- gene.ids
names(q.value_atg8) <- gene.ids
base.mean_atg8 <- res_atg8$baseMean
names(base.mean_atg8) <- gene.ids
# Take a look at upregulated and downregulated genes:
top_20_overexpressed_genes_atg8 <- res_atg8[order(res_atg8$log2FoldChange, decreasing = T),]
kable(top_20_overexpressed_genes_atg8[1:20,-(3:4)])
|
|
baseMean
|
log2FoldChange
|
pvalue
|
padj
|
|
Cre01.g031700.v5.5
|
71.40996
|
22.31599
|
0e+00
|
0.00e+00
|
|
Cre03.g169150.v5.5
|
84.33879
|
22.21561
|
0e+00
|
0.00e+00
|
|
Cre03.g206900.v5.5
|
39.26690
|
22.20814
|
0e+00
|
0.00e+00
|
|
Cre05.g242750.v5.5
|
31.80983
|
22.17803
|
0e+00
|
0.00e+00
|
|
Cre03.g183150.v5.5
|
52.34164
|
22.08262
|
0e+00
|
0.00e+00
|
|
Cre16.g650000.v5.5
|
31.64993
|
21.81398
|
0e+00
|
2.00e-07
|
|
Cre01.g055457.v5.5
|
18.41569
|
21.70069
|
3e-07
|
8.00e-06
|
|
Cre06.g282826.v5.5
|
42.54245
|
21.65897
|
0e+00
|
0.00e+00
|
|
Cre16.g688900.v5.5
|
26.82071
|
21.59666
|
0e+00
|
3.00e-07
|
|
Cre24.g755097.v5.5
|
25.65882
|
21.46373
|
0e+00
|
4.00e-07
|
|
Cre06.g305302.v5.5
|
18.84289
|
21.35537
|
5e-07
|
1.18e-05
|
|
Cre07.g346900.v5.5
|
87.77094
|
21.16046
|
0e+00
|
0.00e+00
|
|
Cre01.g009650.v5.5
|
19.95885
|
21.02765
|
0e+00
|
2.00e-07
|
|
Cre13.g581700.v5.5
|
50.68467
|
20.99081
|
0e+00
|
0.00e+00
|
|
Cre13.g569000.v5.5
|
31.53099
|
20.90911
|
0e+00
|
0.00e+00
|
|
Cre15.g637315.v5.5
|
10.27450
|
20.87907
|
9e-07
|
1.98e-05
|
|
Cre12.g544113.v5.5
|
21.59394
|
20.87294
|
0e+00
|
2.00e-07
|
|
Cre17.g724873.v5.5
|
12.39615
|
20.76123
|
1e-06
|
2.25e-05
|
|
Cre04.g216700.v5.5
|
20.06815
|
20.75520
|
0e+00
|
3.00e-07
|
|
Cre02.g109900.v5.5
|
70.00800
|
20.72860
|
0e+00
|
0.00e+00
|
top_20_overrepressed_genes_atg8 <- res_atg8[order(res_atg8$log2FoldChange, decreasing = F),]
kable(top_20_overrepressed_genes_atg8[1:20,-(3:4)])
|
|
baseMean
|
log2FoldChange
|
pvalue
|
padj
|
|
Cre06.g283000.v5.5
|
20.54063
|
-29.73628
|
0e+00
|
0.0e+00
|
|
Cre17.g735750.v5.5
|
27.85918
|
-29.05778
|
0e+00
|
0.0e+00
|
|
Cre06.g292249.v5.5
|
18.96145
|
-26.98877
|
0e+00
|
0.0e+00
|
|
Cre12.g500250.v5.5
|
28.14109
|
-26.72339
|
0e+00
|
0.0e+00
|
|
Cre02.g115400.v5.5
|
20.74418
|
-25.68656
|
0e+00
|
0.0e+00
|
|
Cre16.g684715.v5.5
|
14.85207
|
-25.34339
|
0e+00
|
1.0e-07
|
|
Cre06.g294250.v5.5
|
123.06077
|
-24.04041
|
0e+00
|
0.0e+00
|
|
Cre08.g376740.v5.5
|
77.07946
|
-23.39426
|
0e+00
|
0.0e+00
|
|
Cre03.g205800.v5.5
|
69.35132
|
-23.21297
|
0e+00
|
1.3e-06
|
|
Cre13.g604550.v5.5
|
56.05740
|
-23.16705
|
0e+00
|
0.0e+00
|
|
Cre07.g315050.v5.5
|
88.92973
|
-23.13839
|
0e+00
|
0.0e+00
|
|
Cre09.g412450.v5.5
|
61.27340
|
-23.03721
|
0e+00
|
0.0e+00
|
|
Cre10.g428000.v5.5
|
52.12164
|
-23.01926
|
0e+00
|
0.0e+00
|
|
Cre13.g592350.v5.5
|
49.26531
|
-22.72924
|
0e+00
|
0.0e+00
|
|
Cre07.g325713.v5.5
|
47.76235
|
-22.60765
|
0e+00
|
0.0e+00
|
|
Cre06.g278168.v5.5
|
57.62603
|
-22.41235
|
0e+00
|
0.0e+00
|
|
Cre12.g547727.v5.5
|
29.91381
|
-22.32585
|
0e+00
|
0.0e+00
|
|
Cre05.g246250.v5.5
|
43.36180
|
-22.26541
|
0e+00
|
0.0e+00
|
|
Cre01.g010832.v5.5
|
37.26769
|
-22.11360
|
0e+00
|
0.0e+00
|
|
Cre16.g695200.v5.5
|
28.69367
|
-22.10441
|
1e-07
|
2.5e-06
|
activated.genes.deseq2_atg8 <- gene.ids[log.fold.change_atg8 > 1 & q.value_atg8 < 0.05]
activated.genes.deseq2_atg8 <- activated.genes.deseq2_atg8[!is.na(activated.genes.deseq2_atg8)]
activated.genes.deseq2_atg8_good <- substr(activated.genes.deseq2_atg8,1,nchar(activated.genes.deseq2_atg8)-5)
write.table(activated.genes.deseq2_atg8_good,file="activated.genes.deseq2_atg8.txt",sep="\t",quote=F,row.names = F)
repressed.genes.deseq2_atg8 <- gene.ids[log.fold.change_atg8 < -1 & q.value_atg8 < 0.05]
repressed.genes.deseq2_atg8 <- repressed.genes.deseq2_atg8[!is.na(repressed.genes.deseq2_atg8)]
repressed.genes.deseq2_atg8_good <- substr(repressed.genes.deseq2_atg8,1,nchar(repressed.genes.deseq2_atg8)-5)
write.table(repressed.genes.deseq2_atg8_good,file="repressed.genes.deseq2_atg8.txt",sep="\t",quote=F,row.names = F)
length(activated.genes.deseq2_atg8_good)
## [1] 618
length(repressed.genes.deseq2_atg8_good)
## [1] 441
MA plot
plotMA(res_atg8,alpha=0.05, ylim=c(-30,30),colSig="grey")
abline(h=c(-1,1), col="black", lwd=1)
points(x = base.meanWT[activated.genes.deseq2_atg8],
y = log.fold.change_atg8[activated.genes.deseq2_atg8],col="red",cex=0.2,pch=19)
points(x = base.mean_atg8[repressed.genes.deseq2_atg8],
y = log.fold.change_atg8[repressed.genes.deseq2_atg8],col="blue",cex=0.2,pch=19)
text(x =100000, y = 15, label = length(activated.genes.deseq2_atg8_good),
col = "red", # Color del texto
font = 2, # Estilo
cex = 1.5) # Tamaño
text(x =100000, y = -15, label = length(repressed.genes.deseq2_atg8_good),
col = "blue", # Color del texto
font = 2, # Estilo
cex = 1.5) # Tamaño
legend("bottomright", legend = c("Overexpressed", "Repressed"),
lwd = 3, col = c("red", "blue"))

Plotting log2 fold change with target genes for wt and atg8 treated with cerulenin
# Test if lists of target genes are in the filtered gene.expression matrix
intersect(ATG_markers$PhyCode, gene.ids) # Check
## [1] "Cre09.g391245.v5.5" "Cre01.g045600.v5.5" "Cre02.g102350.v5.5"
## [4] "Cre12.g510100.v5.5" "Cre14.g630907.v5.5" "Cre03.g165215.v5.5"
## [7] "Cre16.g689650.v5.5" "Cre09.g391500.v5.5" "Cre12.g557000.v5.5"
## [10] "Cre16.g659000.v5.5" "Cre08.g377600.v5.5" "Cre10.g457550.v5.5"
## [13] "Cre07.g334700.v5.5" "Cre06.g290500.v5.5" "Cre01.g035500.v5.5"
intersect(Chloro_stress_markers$PhyCode, gene.ids) # Check
## [1] "Cre13.g583550.v5.5" "Cre11.g468050.v5.5" "Cre06.g250100.v5.5"
## [4] "Cre03.g145787.v5.5" "Cre14.g617450.v5.5" "Cre14.g617400.v5.5"
## [7] "Cre02.g090850.v5.5" "Cre17.g696400.v5.5" "Cre12.g559150.v5.5"
## [10] "Cre07.g341600.v5.5" "Cre12.g507650.v5.5" "Cre01.g009900.v5.5"
## [13] "Cre10.g430400.v5.5" "Cre12.g520600.v5.5"
intersect(Chloro_redox_markers$PhyCode, gene.ids) # Check
## [1] "Cre02.g087700.v5.5" "Cre06.g285150.v5.5" "Cre10.g458450.v5.5"
## [4] "Cre07.g325743.v5.5" "Cre16.g688550.v5.5" "Cre01.g054150.v5.5"
## [7] "Cre10.g422300.v5.5"
ATG Markers
bar.names_ATG <- ATG_markers$Name
bar.phycodes_ATG <- ATG_markers$PhyCode
wt_treatment_ATG <- as.data.frame(resWT) %>% rownames_to_column("GeneID")
wt_treatment_ATG <- wt_treatment_ATG[wt_treatment_ATG$GeneID %in% bar.phycodes_ATG,]
wt_treatment_ATG$GeneID <- ATG_markers$Name[match(wt_treatment_ATG$GeneID,ATG_markers$PhyCode)]
row.names(wt_treatment_ATG) <- wt_treatment_ATG$GeneID
wt_treatment_ATG_fold.change <- wt_treatment_ATG[-c(1:2,4:7)]
colnames(wt_treatment_ATG_fold.change) <- c("WT cerulenin")
wt_treatment_ATG_fold.change.matrix <- as.matrix(wt_treatment_ATG_fold.change)
wt_treatment_ATG$GeneID[wt_treatment_ATG$log2FoldChange >1 & wt_treatment_ATG$padj < 0.05]
## [1] "ATG7" "ATG3" "ATG14" "ATG2" "ATG101" "ATG8"
atg8_treatment_ATG <- as.data.frame(res_atg8) %>% rownames_to_column("GeneID")
atg8_treatment_ATG <- atg8_treatment_ATG[atg8_treatment_ATG$GeneID %in% bar.phycodes_ATG,]
atg8_treatment_ATG$GeneID <- ATG_markers$Name[match(atg8_treatment_ATG$GeneID,ATG_markers$PhyCode)]
row.names(atg8_treatment_ATG) <- atg8_treatment_ATG$GeneID
atg8_treatment_ATG_fold.change <- atg8_treatment_ATG[-c(1:2,4:7)]
colnames(atg8_treatment_ATG_fold.change) <- c("atg8 cerulenin")
atg8_treatment_ATG$GeneID[atg8_treatment_ATG$log2FoldChange >1 & atg8_treatment_ATG$padj < 0.05]
## [1] "VPS34" "ATG7" "ATG3" "ATG14" "ATG2" "ATG101" "ATG8" "ATG13"
atg8_treatment_ATG_fold.change.matrix <- as.matrix(atg8_treatment_ATG_fold.change)
wt_atg8_treatment_ATG_fold.change.matrix <- cbind(wt_treatment_ATG_fold.change.matrix,
atg8_treatment_ATG_fold.change.matrix)
# Transpose matrix:
wt_atg8_treatment_ATG_fold.change.matrix <- t(wt_atg8_treatment_ATG_fold.change.matrix)
# Remove columns where for both genotypes target ATG genes have log2 fold change < 1
# and padj > 0.05
# Colums to remove: VPS15, ATG5, ATG12, ATG9, ATG5, ATG6, ATG1, ATG18, ATG4
wt_atg8_treatment_ATG_fold.change.matrix <- wt_atg8_treatment_ATG_fold.change.matrix[, -c(1,5,7:8,13:16)]
wt_atg8_treatment_ATG_fold.change.matrix<-wt_atg8_treatment_ATG_fold.change.matrix[,order(colnames(wt_atg8_treatment_ATG_fold.change.matrix), decreasing = TRUE)]
barplot(wt_atg8_treatment_ATG_fold.change.matrix, beside = TRUE,
legend.text = TRUE, ylab = "log2 fold change", ylim = c(0,12),
xpd = FALSE, cex.names = 0.5, col=colorRampPalette(rev(brewer.pal(2,"Blues")))(2))
## Warning in brewer.pal(2, "Blues"): minimal value for n is 3, returning requested palette with 3 different levels

Chloroplast chaperones and stress markers
bar.names_ChlStress <- Chloro_stress_markers$Name
bar.phycodes_ChlStress <- Chloro_stress_markers$PhyCode
wt_treatment_ChlStress <- as.data.frame(resWT) %>% rownames_to_column("GeneID")
wt_treatment_ChlStress <- wt_treatment_ChlStress[wt_treatment_ChlStress$GeneID %in% bar.phycodes_ChlStress,]
wt_treatment_ChlStress$GeneID <- Chloro_stress_markers$Name[match(wt_treatment_ChlStress$GeneID,
Chloro_stress_markers$PhyCode)]
row.names(wt_treatment_ChlStress) <- wt_treatment_ChlStress$GeneID
wt_treatment_ChlStress_fold.change <- wt_treatment_ChlStress[-c(1:2,4:7)]
colnames(wt_treatment_ChlStress_fold.change) <- c("WT cerulenin")
wt_treatment_ChlStress_fold.change.matrix <- as.matrix(wt_treatment_ChlStress_fold.change)
wt_treatment_ChlStress$GeneID[wt_treatment_ChlStress$log2FoldChange >1 & wt_treatment_ChlStress$padj < 0.05]
## [1] "VIPP2" "HSP22E" "VIPP1" "HSP22F" "CLPB3"
atg8_treatment_ChlStress <- as.data.frame(res_atg8) %>% rownames_to_column("GeneID")
atg8_treatment_ChlStress <- atg8_treatment_ChlStress[atg8_treatment_ChlStress$GeneID %in% bar.phycodes_ChlStress,]
atg8_treatment_ChlStress$GeneID <- Chloro_stress_markers$Name[match(atg8_treatment_ChlStress$GeneID,
Chloro_stress_markers$PhyCode)]
row.names(atg8_treatment_ChlStress) <- atg8_treatment_ChlStress$GeneID
atg8_treatment_ChlStress_fold.change <- atg8_treatment_ChlStress[-c(1:2,4:7)]
colnames(atg8_treatment_ChlStress_fold.change) <- c("atg8 cerulenin")
atg8_treatment_ChlStress$GeneID[atg8_treatment_ChlStress$log2FoldChange >1 & atg8_treatment_ChlStress$padj < 0.05]
## [1] "VIPP2" "HSP22E" "VIPP1" "HSP22F" "CLPB3"
atg8_treatment_ChlStress_fold.change.matrix <- as.matrix(atg8_treatment_ChlStress_fold.change)
wt_atg8_treatment_ChlStress_fold.change.matrix <- cbind(wt_treatment_ChlStress_fold.change.matrix,
atg8_treatment_ChlStress_fold.change.matrix)
# Transpose matrix:
wt_atg8_treatment_ChlStress_fold.change.matrix <- t(wt_atg8_treatment_ChlStress_fold.change.matrix)
# Remove columns where for both genotypes target ATG genes have log2 fold change < 1
# and padj > 0.05
# Colums to remove: HSP22C, RPL37, CDJ2, CDJ1, CGE1, HSP70B, CDJ3, CLPS1, PRPS6, CLPS2
wt_atg8_treatment_ChlStress_fold.change.matrix <- wt_atg8_treatment_ChlStress_fold.change.matrix[, -c(2:3, 5:10,12,14)]
wt_atg8_treatment_ChlStress_fold.change.matrix<-wt_atg8_treatment_ChlStress_fold.change.matrix[,order(colnames(wt_atg8_treatment_ChlStress_fold.change.matrix), decreasing = TRUE)]
barplot(wt_atg8_treatment_ChlStress_fold.change.matrix, beside = TRUE,
legend.text = TRUE, ylab = "log2 fold change", ylim = c(0,6),
xpd = FALSE, cex.names = 0.5, col=colorRampPalette(rev(brewer.pal(2,"Blues")))(2))
## Warning in brewer.pal(2, "Blues"): minimal value for n is 3, returning requested palette with 3 different levels

Chloroplast redox markers
bar.names_ChlRed <- Chloro_redox_markers$Name
bar.phycodes_ChlRed <- Chloro_redox_markers$PhyCode
wt_treatment_ChlRed <- as.data.frame(resWT) %>% rownames_to_column("GeneID")
wt_treatment_ChlRed <- wt_treatment_ChlRed[wt_treatment_ChlRed$GeneID %in% bar.phycodes_ChlRed,]
wt_treatment_ChlRed$GeneID <- Chloro_redox_markers$Name[match(wt_treatment_ChlRed$GeneID,
Chloro_redox_markers$PhyCode)]
row.names(wt_treatment_ChlRed) <- wt_treatment_ChlRed$GeneID
wt_treatment_ChlRed_fold.change <- wt_treatment_ChlRed[-c(1:2,4:7)]
colnames(wt_treatment_ChlRed_fold.change) <- c("WT cerulenin")
wt_treatment_ChlRed_fold.change.matrix <- as.matrix(wt_treatment_ChlRed_fold.change)
wt_treatment_ChlRed$GeneID[wt_treatment_ChlRed$log2FoldChange >1 & wt_treatment_ChlRed$padj < 0.05]
## [1] "GSTS1" "APX1" "GPX5"
atg8_treatment_ChlRed <- as.data.frame(res_atg8) %>% rownames_to_column("GeneID")
atg8_treatment_ChlRed <- atg8_treatment_ChlRed[atg8_treatment_ChlRed$GeneID %in% bar.phycodes_ChlRed,]
atg8_treatment_ChlRed$GeneID <- Chloro_redox_markers$Name[match(atg8_treatment_ChlRed$GeneID,
Chloro_redox_markers$PhyCode)]
row.names(atg8_treatment_ChlRed) <- atg8_treatment_ChlRed$GeneID
atg8_treatment_ChlRed_fold.change <- atg8_treatment_ChlRed[-c(1:2,4:7)]
colnames(atg8_treatment_ChlRed_fold.change) <- c("atg8 cerulenin")
atg8_treatment_ChlRed$GeneID[atg8_treatment_ChlRed$log2FoldChange >1 & atg8_treatment_ChlRed$padj < 0.05]
## [1] "GSTS1" "APX1" "GPX5"
atg8_treatment_ChlRed_fold.change.matrix <- as.matrix(atg8_treatment_ChlRed_fold.change)
wt_atg8_treatment_ChlRed_fold.change.matrix <- cbind(wt_treatment_ChlRed_fold.change.matrix,
atg8_treatment_ChlRed_fold.change.matrix)
# Transpose matrix:
wt_atg8_treatment_ChlRed_fold.change.matrix <- t(wt_atg8_treatment_ChlRed_fold.change.matrix)
# Remove columns where for both genotypes target ATG genes have log2 fold change < 1
# and padj > 0.05
# Colums to remove: APX2, GRX3, PRX6, NTRC1
wt_atg8_treatment_ChlRed_fold.change.matrix <- wt_atg8_treatment_ChlRed_fold.change.matrix[, -c(3:6)]
wt_atg8_treatment_ChlRed_fold.change.matrix<-wt_atg8_treatment_ChlRed_fold.change.matrix[,order(colnames(wt_atg8_treatment_ChlRed_fold.change.matrix), decreasing = TRUE)]
barplot(wt_atg8_treatment_ChlRed_fold.change.matrix, beside = TRUE,
legend.text = TRUE, ylab = "log2 fold change", ylim = c(0,6),
xpd = FALSE, cex.names = 0.5, col=colorRampPalette(rev(brewer.pal(2,"Blues")))(2))
## Warning in brewer.pal(2, "Blues"): minimal value for n is 3, returning requested palette with 3 different levels

Difference between wild-type and mutant 0h
Gene obtention
res_WT_atg8_0h = results(dds, contrast=c("genotype","atg8","wt"),
cooksCutoff = FALSE, independentFiltering = FALSE)
res_WT_atg8_0h
## log2 fold change (MLE): genotype atg8 vs wt
## Wald test p-value: genotype atg8 vs wt
## DataFrame with 13977 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue
## <numeric> <numeric> <numeric> <numeric> <numeric>
## Cre14.g619950.v5.5 307.8280 -1.0486329 0.719601 -1.457241 0.145050
## Cre01.g048900.v5.5 76.4131 0.2567535 1.830265 0.140282 0.888437
## Cre12.g543400.v5.5 428.9410 0.4253013 0.320785 1.325814 0.184901
## Cre16.g685837.v5.5 255.7288 0.0698002 0.456223 0.152996 0.878402
## Cre06.g278159.v5.5 9770.0557 -0.0796437 0.124401 -0.640217 0.522031
## ... ... ... ... ... ...
## Cre06.g278104.v5.5 889.54623 -0.604567 0.310855 -1.944851 0.0517930
## Cre17.g717150.v5.5 146.81853 0.994081 1.170782 0.849074 0.3958403
## Cre10.g443600.v5.5 170.10265 0.127679 1.079078 0.118322 0.9058126
## Cre14.g621400.v5.5 5.13639 6.786913 4.117852 1.648168 0.0993182
## Cre10.g441750.v5.5 1599.45474 0.269795 0.181866 1.483486 0.1379453
## padj
## <numeric>
## Cre14.g619950.v5.5 0.728842
## Cre01.g048900.v5.5 0.999313
## Cre12.g543400.v5.5 0.799370
## Cre16.g685837.v5.5 0.999313
## Cre06.g278159.v5.5 0.999313
## ... ...
## Cre06.g278104.v5.5 0.445209
## Cre17.g717150.v5.5 0.981999
## Cre10.g443600.v5.5 0.999313
## Cre14.g621400.v5.5 0.619202
## Cre10.g441750.v5.5 0.711462
log.fold.change_WT_atg8_0h <- res_WT_atg8_0h$log2FoldChange
q.value_WT_atg8_0h <- res_WT_atg8_0h$padj
names(log.fold.change_WT_atg8_0h) <- gene.ids
names(q.value_WT_atg8_0h) <- gene.ids
base.mean_WT_atg8_0h <- res_WT_atg8_0h$baseMean
names(base.mean_WT_atg8_0h) <- gene.ids
# Take a look at upregulated and downregulated genes:
top_20_overexpressed_genes_WT_atg8_0h <- res_WT_atg8_0h[order(res_WT_atg8_0h$log2FoldChange, decreasing = T),]
kable(top_20_overexpressed_genes_WT_atg8_0h[1:20,-(3:4)])
|
|
baseMean
|
log2FoldChange
|
pvalue
|
padj
|
|
Cre10.g466200.v5.5
|
15.235117
|
20.585945
|
0.0000013
|
0.0001480
|
|
Cre07.g331350.v5.5
|
19.178658
|
19.940809
|
0.0000005
|
0.0000721
|
|
Cre02.g143107.v5.5
|
7.114070
|
19.404147
|
0.0000051
|
0.0004478
|
|
Cre10.g465763.v5.5
|
17.445204
|
18.935435
|
0.0000086
|
0.0006709
|
|
Cre06.g250650.v5.5
|
13.336058
|
14.980359
|
0.0004177
|
0.0159081
|
|
Cre13.g567450.v5.5
|
47.816726
|
9.974608
|
0.0028980
|
0.0675091
|
|
Cre06.g310900.v5.5
|
12.697053
|
8.158493
|
0.0541851
|
0.4562322
|
|
Cre03.g167622.v5.5
|
13.099835
|
8.019118
|
0.0284711
|
0.3099224
|
|
Cre08.g358565.v5.5
|
7.390831
|
7.384139
|
0.0824190
|
0.5624718
|
|
Cre02.g120200.v5.5
|
7.613366
|
7.380495
|
0.0815050
|
0.5587032
|
|
Cre15.g635100.v5.5
|
16.735402
|
7.208932
|
0.0296900
|
0.3160524
|
|
Cre17.g713500.v5.5
|
5.698800
|
7.072735
|
0.0671212
|
0.5118263
|
|
Cre13.g569900.v5.5
|
12.543505
|
6.987735
|
0.1003772
|
0.6221609
|
|
Cre09.g404552.v5.5
|
7.466950
|
6.930540
|
0.1032099
|
0.6306568
|
|
Cre14.g621400.v5.5
|
5.136393
|
6.786913
|
0.0993182
|
0.6192016
|
|
Cre11.g480502.v5.5
|
4.849025
|
6.781913
|
0.0473071
|
0.4194994
|
|
Cre09.g395925.v5.5
|
10.571970
|
6.551579
|
0.1236302
|
0.6819178
|
|
Cre03.g144444.v5.5
|
8.399966
|
6.412211
|
0.1318948
|
0.6990198
|
|
Cre06.g265200.v5.5
|
3.559183
|
6.353395
|
0.1041611
|
0.6328170
|
|
Cre02.g091226.v5.5
|
16.082502
|
6.198403
|
0.1454064
|
0.7288417
|
top_20_overrepressed_genes_WT_atg8_0h <- res_WT_atg8_0h[order(res_WT_atg8_0h$log2FoldChange, decreasing = F),]
kable(top_20_overrepressed_genes_WT_atg8_0h[1:20,-(3:4)])
|
|
baseMean
|
log2FoldChange
|
pvalue
|
padj
|
|
Cre07.g346900.v5.5
|
87.770943
|
-23.11461
|
0.0e+00
|
0.0000000
|
|
Cre08.g364650.v5.5
|
58.216483
|
-22.48245
|
0.0e+00
|
0.0000000
|
|
Cre02.g109900.v5.5
|
70.008000
|
-22.42246
|
0.0e+00
|
0.0000000
|
|
Cre03.g183150.v5.5
|
52.341644
|
-21.67106
|
0.0e+00
|
0.0000000
|
|
Cre10.g433650.v5.5
|
19.955260
|
-21.51465
|
1.0e-07
|
0.0000204
|
|
Cre02.g077951.v5.5
|
22.406577
|
-21.47873
|
1.0e-07
|
0.0000134
|
|
Cre03.g169150.v5.5
|
84.338786
|
-21.16178
|
0.0e+00
|
0.0000000
|
|
Cre02.g095141.v5.5
|
19.617552
|
-21.06563
|
3.0e-07
|
0.0000407
|
|
Cre13.g581700.v5.5
|
50.684675
|
-21.00907
|
0.0e+00
|
0.0000000
|
|
Cre10.g452500.v5.5
|
27.671211
|
-20.91627
|
5.0e-07
|
0.0000661
|
|
Cre06.g282826.v5.5
|
42.542451
|
-20.82930
|
0.0e+00
|
0.0000000
|
|
Cre12.g494600.v5.5
|
20.255885
|
-20.54698
|
3.0e-07
|
0.0000407
|
|
Cre07.g350050.v5.5
|
12.690641
|
-20.50715
|
1.4e-06
|
0.0001607
|
|
Cre01.g031700.v5.5
|
71.409961
|
-20.47323
|
0.0e+00
|
0.0000000
|
|
Cre03.g206900.v5.5
|
39.266900
|
-20.37765
|
0.0e+00
|
0.0000001
|
|
Cre09.g395732.v5.5
|
26.133934
|
-20.36151
|
1.0e-07
|
0.0000217
|
|
Cre06.g297049.v5.5
|
9.569606
|
-20.31073
|
1.2e-06
|
0.0001435
|
|
Cre17.g698950.v5.5
|
9.775073
|
-20.25779
|
1.9e-06
|
0.0002031
|
|
Cre03.g191650.v5.5
|
21.334812
|
-20.15767
|
0.0e+00
|
0.0000030
|
|
Cre16.g669350.v5.5
|
14.906514
|
-20.03924
|
9.0e-07
|
0.0001214
|
activated.genes.deseq2_WT_atg8_0h <- gene.ids[log.fold.change_WT_atg8_0h > 1 & q.value_WT_atg8_0h < 0.05]
activated.genes.deseq2_WT_atg8_0h <- activated.genes.deseq2_WT_atg8_0h[!is.na(activated.genes.deseq2_WT_atg8_0h)]
activated.genes.deseq2_WT_atg8_0h_good <- substr(activated.genes.deseq2_WT_atg8_0h,1,nchar(activated.genes.deseq2_WT_atg8_0h)-5)
write.table(activated.genes.deseq2_WT_atg8_0h_good,file="activated.genes.deseq2_WT_atg8_0h.txt",sep="\t",quote=F,row.names = F)
repressed.genes.deseq2_WT_atg8_0h <- gene.ids[log.fold.change_WT_atg8_0h < -1 & q.value_WT_atg8_0h < 0.05]
repressed.genes.deseq2_WT_atg8_0h <- repressed.genes.deseq2_WT_atg8_0h[!is.na(repressed.genes.deseq2_WT_atg8_0h)]
repressed.genes.deseq2_WT_atg8_0h_good <- substr(repressed.genes.deseq2_WT_atg8_0h,1,nchar(repressed.genes.deseq2_WT_atg8_0h)-5)
write.table(repressed.genes.deseq2_WT_atg8_0h_good,file="repressed.genes.deseq2_WT_atg8_0h.txt",sep="\t",quote=F,row.names = F)
length(activated.genes.deseq2_WT_atg8_0h_good)
## [1] 83
length(repressed.genes.deseq2_WT_atg8_0h_good)
## [1] 117
MA plot
plotMA(res_WT_atg8_0h,alpha=0.05, ylim=c(-30,30),colSig="grey")
abline(h=c(-1,1), col="black", lwd=1)
points(x = base.mean_WT_atg8_0h[activated.genes.deseq2_WT_atg8_0h],
y = log.fold.change_WT_atg8_0h[activated.genes.deseq2_WT_atg8_0h],col="red",cex=0.2,pch=19)
points(x = base.mean_WT_atg8_0h[repressed.genes.deseq2_WT_atg8_0h],
y = log.fold.change_WT_atg8_0h[repressed.genes.deseq2_WT_atg8_0h],col="blue",cex=0.2,pch=19)
text(x =100000, y = 15, label = length(activated.genes.deseq2_WT_atg8_0h_good),
col = "red", # Color del texto
font = 2, # Estilo
cex = 1.5) # Tamaño
text(x =100000, y = -15, label = length(repressed.genes.deseq2_WT_atg8_0h_good),
col = "blue", # Color del texto
font = 2, # Estilo
cex = 1.5) # Tamaño
legend("bottomright", legend = c("Overexpressed", "Repressed"),
lwd = 3, col = c("red", "blue"))

Difference between wild-type and mutant 4h
Gene obtention
res_WT_atg8_4h = results(dds, list( c("genotype_atg8_vs_wt","genotypeatg8.condition4h")),
cooksCutoff = FALSE, independentFiltering = FALSE )
res_WT_atg8_4h
## log2 fold change (MLE): genotype_atg8_vs_wt+genotypeatg8.condition4h effect
## Wald test p-value: genotype_atg8_vs_wt+genotypeatg8.condition4h effect
## DataFrame with 13977 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue
## <numeric> <numeric> <numeric> <numeric> <numeric>
## Cre14.g619950.v5.5 307.8280 -1.094574 0.719128 -1.522086 0.12798748
## Cre01.g048900.v5.5 76.4131 -0.476455 1.829920 -0.260369 0.79457914
## Cre12.g543400.v5.5 428.9410 0.631098 0.320734 1.967667 0.04910641
## Cre16.g685837.v5.5 255.7288 0.130937 0.454228 0.288264 0.77314498
## Cre06.g278159.v5.5 9770.0557 0.361482 0.124145 2.911778 0.00359378
## ... ... ... ... ... ...
## Cre06.g278104.v5.5 889.54623 0.228314 0.310320 0.735738 0.4618902
## Cre17.g717150.v5.5 146.81853 -0.644081 1.173279 -0.548958 0.5830341
## Cre10.g443600.v5.5 170.10265 1.887998 1.082569 1.743998 0.0811594
## Cre14.g621400.v5.5 5.13639 4.605473 4.141349 1.112071 0.2661077
## Cre10.g441750.v5.5 1599.45474 0.277861 0.181312 1.532499 0.1253994
## padj
## <numeric>
## Cre14.g619950.v5.5 0.5198724
## Cre01.g048900.v5.5 0.9823882
## Cre12.g543400.v5.5 0.3269940
## Cre16.g685837.v5.5 0.9793068
## Cre06.g278159.v5.5 0.0577359
## ... ...
## Cre06.g278104.v5.5 0.879075
## Cre17.g717150.v5.5 0.930899
## Cre10.g443600.v5.5 0.420310
## Cre14.g621400.v5.5 0.733029
## Cre10.g441750.v5.5 0.515712
log.fold.change_WT_atg8_4h <- res_WT_atg8_4h$log2FoldChange
q.value_WT_atg8_4h <- res_WT_atg8_4h$padj
names(log.fold.change_WT_atg8_4h) <- gene.ids
names(q.value_WT_atg8_4h) <- gene.ids
base.mean_WT_atg8_4h <- res_WT_atg8_4h$baseMean
names(base.mean_WT_atg8_4h) <- gene.ids
# Take a look at upregulated and downregulated genes:
top_20_overexpressed_genes_WT_atg8_4h <- res_WT_atg8_4h[order(res_WT_atg8_4h$log2FoldChange, decreasing = T),]
kable(top_20_overexpressed_genes_WT_atg8_4h[1:20,-(3:4)])
|
|
baseMean
|
log2FoldChange
|
pvalue
|
padj
|
|
Cre17.g712700.v5.5
|
42.409362
|
23.20189
|
0e+00
|
0.00e+00
|
|
Cre16.g658400.v5.5
|
59.759725
|
23.01070
|
0e+00
|
0.00e+00
|
|
Cre12.g524450.v5.5
|
32.024953
|
21.89046
|
0e+00
|
4.00e-07
|
|
Cre13.g588250.v5.5
|
49.571334
|
21.80803
|
0e+00
|
0.00e+00
|
|
Cre16.g687966.v5.5
|
33.784426
|
21.79817
|
0e+00
|
0.00e+00
|
|
Cre03.g143827.v5.5
|
66.501562
|
21.66314
|
0e+00
|
0.00e+00
|
|
Cre15.g636600.v5.5
|
58.721594
|
21.62000
|
0e+00
|
0.00e+00
|
|
Cre07.g336500.v5.5
|
36.496724
|
21.41570
|
0e+00
|
4.00e-07
|
|
Cre01.g043250.v5.5
|
102.227854
|
21.36041
|
0e+00
|
0.00e+00
|
|
Cre03.g195410.v5.5
|
16.175576
|
21.30459
|
0e+00
|
5.00e-07
|
|
Cre12.g524350.v5.5
|
18.072167
|
21.27414
|
6e-07
|
3.66e-05
|
|
Cre03.g179400.v5.5
|
39.299172
|
21.24146
|
0e+00
|
0.00e+00
|
|
Cre11.g467532.v5.5
|
33.222102
|
21.23809
|
0e+00
|
0.00e+00
|
|
Cre05.g240100.v5.5
|
25.335507
|
21.21101
|
0e+00
|
3.60e-06
|
|
Cre05.g239000.v5.5
|
46.353687
|
21.15574
|
0e+00
|
0.00e+00
|
|
Cre06.g292600.v5.5
|
15.273678
|
21.14644
|
7e-07
|
4.17e-05
|
|
Cre01.g034380.v5.5
|
38.571871
|
21.10793
|
0e+00
|
0.00e+00
|
|
Cre10.g422350.v5.5
|
8.728187
|
21.10641
|
7e-07
|
4.32e-05
|
|
Cre09.g402050.v5.5
|
54.881325
|
21.10392
|
0e+00
|
0.00e+00
|
|
Cre14.g629840.v5.5
|
37.019186
|
21.03051
|
0e+00
|
0.00e+00
|
top_20_overrepressed_genes_WT_atg8_4h <- res_WT_atg8_4h[order(res_WT_atg8_4h$log2FoldChange, decreasing = F),]
kable(top_20_overrepressed_genes_WT_atg8_4h[1:20,-(3:4)])
|
|
baseMean
|
log2FoldChange
|
pvalue
|
padj
|
|
Cre02.g077951.v5.5
|
22.40658
|
-27.39922
|
0e+00
|
0.00e+00
|
|
Cre10.g433650.v5.5
|
19.95526
|
-26.86531
|
0e+00
|
0.00e+00
|
|
Cre07.g350050.v5.5
|
12.69064
|
-26.22043
|
0e+00
|
1.00e-07
|
|
Cre02.g095141.v5.5
|
19.61755
|
-26.18287
|
0e+00
|
0.00e+00
|
|
Cre07.g315050.v5.5
|
88.92973
|
-23.62099
|
0e+00
|
0.00e+00
|
|
Cre16.g684861.v5.5
|
36.98465
|
-22.61599
|
0e+00
|
0.00e+00
|
|
Cre07.g327079.v5.5
|
23.64919
|
-22.52874
|
1e-07
|
6.40e-06
|
|
Cre13.g592350.v5.5
|
49.26531
|
-22.51609
|
0e+00
|
0.00e+00
|
|
Cre14.g623300.v5.5
|
27.93290
|
-22.39287
|
0e+00
|
0.00e+00
|
|
Cre13.g604550.v5.5
|
56.05740
|
-22.25126
|
0e+00
|
0.00e+00
|
|
Cre02.g118450.v5.5
|
30.20969
|
-22.24960
|
0e+00
|
4.50e-06
|
|
Cre03.g205800.v5.5
|
69.35132
|
-22.19388
|
2e-07
|
1.33e-05
|
|
Cre12.g509600.v5.5
|
38.66322
|
-22.15181
|
0e+00
|
0.00e+00
|
|
Cre06.g278168.v5.5
|
57.62603
|
-22.07722
|
0e+00
|
0.00e+00
|
|
Cre03.g153250.v5.5
|
28.01849
|
-22.07050
|
1e-07
|
9.40e-06
|
|
Cre07.g314700.v5.5
|
37.60714
|
-21.95988
|
0e+00
|
0.00e+00
|
|
Cre16.g678997.v5.5
|
30.44633
|
-21.93106
|
0e+00
|
0.00e+00
|
|
Cre06.g294250.v5.5
|
123.06077
|
-21.90895
|
0e+00
|
0.00e+00
|
|
Cre15.g636250.v5.5
|
24.67298
|
-21.89186
|
0e+00
|
1.80e-06
|
|
Cre01.g010832.v5.5
|
37.26769
|
-21.87643
|
0e+00
|
0.00e+00
|
activated.genes.deseq2_WT_atg8_4h <- gene.ids[log.fold.change_WT_atg8_4h > 1 & q.value_WT_atg8_4h < 0.05]
activated.genes.deseq2_WT_atg8_4h <- activated.genes.deseq2_WT_atg8_4h[!is.na(activated.genes.deseq2_WT_atg8_4h)]
activated.genes.deseq2_WT_atg8_4h_good <- substr(activated.genes.deseq2_WT_atg8_4h,1,nchar(activated.genes.deseq2_WT_atg8_4h)-5)
write.table(activated.genes.deseq2_WT_atg8_4h_good,file="activated.genes.deseq2_WT_atg8_4h.txt",sep="\t",quote=F,row.names = F)
repressed.genes.deseq2_WT_atg8_4h <- gene.ids[log.fold.change_WT_atg8_4h < -1 & q.value_WT_atg8_4h < 0.05]
repressed.genes.deseq2_WT_atg8_4h <- repressed.genes.deseq2_WT_atg8_4h[!is.na(repressed.genes.deseq2_WT_atg8_4h)]
repressed.genes.deseq2_WT_atg8_4h_good <- substr(repressed.genes.deseq2_WT_atg8_4h,1,nchar(repressed.genes.deseq2_WT_atg8_4h)-5)
write.table(repressed.genes.deseq2_WT_atg8_4h_good,file="repressed.genes.deseq2_WT_atg8_4h.txt",sep="\t",quote=F,row.names = F)
length(activated.genes.deseq2_WT_atg8_4h)
## [1] 193
length(repressed.genes.deseq2_WT_atg8_4h)
## [1] 187
MA plot
plotMA(res_WT_atg8_4h, ylim=c(-30,30),colSig="grey")
abline(h=c(-1,1), col="black", lwd=1)
points(x = base.mean_WT_atg8_4h[activated.genes.deseq2_WT_atg8_4h],
y = log.fold.change_WT_atg8_4h[activated.genes.deseq2_WT_atg8_4h],col="red",cex=0.2,pch=19)
points(x = base.mean_WT_atg8_4h[repressed.genes.deseq2_WT_atg8_4h],
y = log.fold.change_WT_atg8_4h[repressed.genes.deseq2_WT_atg8_4h],col="blue",cex=0.2,pch=19)
text(x =100000, y = 15, label = length(activated.genes.deseq2_WT_atg8_4h_good),
col = "red", # Color del texto
font = 2, # Estilo
cex = 1.5) # Tamaño
text(x =100000, y = -15, label = length(repressed.genes.deseq2_WT_atg8_4h_good),
col = "blue", # Color del texto
font = 2, # Estilo
cex = 1.5) # Tamaño
legend("bottomright", legend = c("Overexpressed", "Repressed"),
lwd = 3, col = c("red", "blue"))

Different response for both genotypes (the interaction term)
Gene obtention
res_interaction = results(dds, name="genotypeatg8.condition4h",
cooksCutoff = FALSE, independentFiltering = FALSE)
res_interaction
## log2 fold change (MLE): genotypeatg8.condition4h
## Wald test p-value: genotypeatg8.condition4h
## DataFrame with 13977 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue
## <numeric> <numeric> <numeric> <numeric> <numeric>
## Cre14.g619950.v5.5 307.8280 -0.0459415 1.017335 -0.0451587 0.9639808
## Cre01.g048900.v5.5 76.4131 -0.7332081 2.588131 -0.2832963 0.7769497
## Cre12.g543400.v5.5 428.9410 0.2057964 0.453622 0.4536733 0.6500640
## Cre16.g685837.v5.5 255.7288 0.0611372 0.643788 0.0949648 0.9243428
## Cre06.g278159.v5.5 9770.0557 0.4411260 0.175749 2.5099837 0.0120737
## ... ... ... ... ... ...
## Cre06.g278104.v5.5 889.54623 0.83288039 0.439237 1.896199 0.0579338
## Cre17.g717150.v5.5 146.81853 -1.63816151 1.657499 -0.988333 0.3229896
## Cre10.g443600.v5.5 170.10265 1.76031972 1.528515 1.151653 0.2494635
## Cre14.g621400.v5.5 5.13639 -2.18143969 5.840032 -0.373532 0.7087524
## Cre10.g441750.v5.5 1599.45474 0.00806577 0.256806 0.031408 0.9749441
## padj
## <numeric>
## Cre14.g619950.v5.5 0.999419
## Cre01.g048900.v5.5 0.999419
## Cre12.g543400.v5.5 0.999419
## Cre16.g685837.v5.5 0.999419
## Cre06.g278159.v5.5 0.366382
## ... ...
## Cre06.g278104.v5.5 0.772653
## Cre17.g717150.v5.5 0.999419
## Cre10.g443600.v5.5 0.999419
## Cre14.g621400.v5.5 0.999419
## Cre10.g441750.v5.5 0.999419
log.fold.change_interaction <- res_interaction$log2FoldChange
q.value_interaction <- res_interaction$padj
names(log.fold.change_interaction) <- gene.ids
names(q.value_interaction) <- gene.ids
base.mean_interaction <- res_interaction$baseMean
names(base.mean_interaction) <- gene.ids
# Take a look at upregulated and downregulated genes:
top_20_overexpressed_genes_interaction <- res_interaction[order(res_interaction$log2FoldChange, decreasing = T),]
kable(top_20_overexpressed_genes_interaction[1:20,-(3:4)])
|
|
baseMean
|
log2FoldChange
|
pvalue
|
padj
|
|
Cre04.g214769.v5.5
|
8.578910
|
30.00000
|
0.0000006
|
0.0001753
|
|
Cre07.g331114.v5.5
|
14.030683
|
30.00000
|
0.0000006
|
0.0001753
|
|
Cre03.g196000.v5.5
|
10.329530
|
30.00000
|
0.0000006
|
0.0001753
|
|
Cre04.g213950.v5.5
|
7.990901
|
30.00000
|
0.0000006
|
0.0001753
|
|
Cre06.g292450.v5.5
|
8.428540
|
30.00000
|
0.0000006
|
0.0001753
|
|
Cre11.g481104.v5.5
|
17.175085
|
30.00000
|
0.0000003
|
0.0001134
|
|
Cre04.g219750.v5.5
|
10.123388
|
30.00000
|
0.0000006
|
0.0001753
|
|
Cre02.g087500.v5.5
|
10.921086
|
30.00000
|
0.0000006
|
0.0001753
|
|
Cre14.g633900.v5.5
|
24.339891
|
30.00000
|
0.0000001
|
0.0000342
|
|
Cre02.g099100.v5.5
|
33.227640
|
30.00000
|
0.0000000
|
0.0000036
|
|
Cre16.g658300.v5.5
|
46.217102
|
23.52922
|
0.0000000
|
0.0000006
|
|
Cre03.g183150.v5.5
|
52.341644
|
22.83219
|
0.0000000
|
0.0000245
|
|
Cre03.g206900.v5.5
|
39.266900
|
22.28228
|
0.0000005
|
0.0001753
|
|
Cre07.g319600.v5.5
|
9.767238
|
22.09547
|
0.0001150
|
0.0143253
|
|
Cre02.g098550.v5.5
|
10.725460
|
22.02765
|
0.0001957
|
0.0199689
|
|
Cre17.g699550.v5.5
|
11.544327
|
22.01875
|
0.0001931
|
0.0199689
|
|
Cre07.g346900.v5.5
|
87.770943
|
21.99156
|
0.0000000
|
0.0000000
|
|
Cre07.g330150.v5.5
|
6.660717
|
21.94332
|
0.0002175
|
0.0218705
|
|
Cre03.g179400.v5.5
|
39.299172
|
21.70125
|
0.0000006
|
0.0001753
|
|
Cre16.g658400.v5.5
|
59.759725
|
21.68779
|
0.0000001
|
0.0000285
|
top_20_overrepressed_genes_interaction <- res_interaction[order(res_interaction$log2FoldChange, decreasing = F),]
kable(top_20_overrepressed_genes_interaction[1:20,-(3:4)])
|
|
baseMean
|
log2FoldChange
|
pvalue
|
padj
|
|
Cre06.g250650.v5.5
|
13.336058
|
-30.00000
|
0.0000006
|
0.0001753
|
|
Cre16.g691351.v5.5
|
18.498013
|
-24.96833
|
0.0000052
|
0.0012079
|
|
Cre07.g315050.v5.5
|
88.929727
|
-24.78043
|
0.0000000
|
0.0000000
|
|
Cre02.g118450.v5.5
|
30.209685
|
-24.24031
|
0.0000181
|
0.0034299
|
|
Cre13.g592350.v5.5
|
49.265306
|
-24.11483
|
0.0000000
|
0.0000053
|
|
Cre13.g604550.v5.5
|
56.057399
|
-23.90740
|
0.0000000
|
0.0000001
|
|
Cre16.g678997.v5.5
|
30.446334
|
-23.88634
|
0.0000003
|
0.0001089
|
|
Cre16.g695200.v5.5
|
28.693675
|
-23.54711
|
0.0000408
|
0.0066344
|
|
Cre12.g544327.v5.5
|
9.507146
|
-23.34612
|
0.0000802
|
0.0112272
|
|
Cre09.g413000.v5.5
|
10.374182
|
-23.30969
|
0.0000816
|
0.0112272
|
|
Cre03.g153250.v5.5
|
28.018485
|
-23.29247
|
0.0000548
|
0.0086987
|
|
Cre16.g684861.v5.5
|
36.984651
|
-23.08733
|
0.0000000
|
0.0000077
|
|
Cre12.g486450.v5.5
|
11.227031
|
-22.92555
|
0.0001045
|
0.0133968
|
|
Cre08.g382825.v5.5
|
12.828425
|
-22.85710
|
0.0001092
|
0.0138783
|
|
Cre14.g609150.v5.5
|
11.132919
|
-22.63246
|
0.0001300
|
0.0152511
|
|
Cre03.g151300.v5.5
|
13.711523
|
-22.62959
|
0.0001261
|
0.0150679
|
|
Cre03.g205800.v5.5
|
69.351320
|
-22.56410
|
0.0001282
|
0.0151822
|
|
Cre10.g466200.v5.5
|
15.235117
|
-22.51100
|
0.0001425
|
0.0161548
|
|
Cre17.g736350.v5.5
|
15.843663
|
-22.49394
|
0.0001378
|
0.0159119
|
|
Cre04.g215001.v5.5
|
11.408657
|
-22.48296
|
0.0001406
|
0.0161085
|
activated.genes.deseq2_interaction <- gene.ids[log.fold.change_interaction > 1 & q.value_interaction < 0.05]
activated.genes.deseq2_interaction <- activated.genes.deseq2_interaction[!is.na(activated.genes.deseq2_interaction)]
activated.genes.deseq2_interaction_good <- substr(activated.genes.deseq2_interaction,1,nchar(activated.genes.deseq2_interaction)-5)
write.table(activated.genes.deseq2_interaction_good,file="activated.genes.deseq2_interaction.txt",sep="\t",quote=F,row.names = F)
repressed.genes.deseq2_interaction <- gene.ids[log.fold.change_interaction < -1 & q.value_interaction < 0.05]
repressed.genes.deseq2_interaction <- repressed.genes.deseq2_interaction[!is.na(repressed.genes.deseq2_interaction)]
repressed.genes.deseq2_interaction_good <- substr(repressed.genes.deseq2_interaction,1,nchar(repressed.genes.deseq2_interaction)-5)
write.table(repressed.genes.deseq2_interaction_good,file="repressed.genes.deseq2_interaction.txt",sep="\t",quote=F,row.names = F)
length(activated.genes.deseq2_interaction)
## [1] 103
length(repressed.genes.deseq2_interaction)
## [1] 87
MA plot
plotMA(res_interaction, ylim=c(-30,30),colSig="grey")
abline(h=c(-1,1), col="black", lwd=1)
points(x = base.mean_interaction[activated.genes.deseq2_interaction],
y = log.fold.change_interaction[activated.genes.deseq2_interaction],col="red",cex=0.2,pch=19)
points(x = base.mean_interaction[repressed.genes.deseq2_interaction],
y = log.fold.change_interaction[repressed.genes.deseq2_interaction],col="blue",cex=0.2,pch=19)
text(x =100000, y = 15, label = length(activated.genes.deseq2_interaction_good),
col = "red", # Color del texto
font = 2, # Estilo
cex = 1.5) # Tamaño
text(x =100000, y = -15, label = length(repressed.genes.deseq2_interaction_good),
col = "blue", # Color del texto
font = 2, # Estilo
cex = 1.5) # Tamaño
legend("bottomright", legend = c("Overexpressed", "Repressed"),
lwd = 3, col = c("red", "blue"))

Heatmap
interaction_sigs <- res_interaction[res_interaction$padj < 0.05,]
interaction_df <- as.data.frame(interaction_sigs)
interaction_dftop <- interaction_df[interaction_df$baseMean > 50 &
abs(interaction_df$log2FoldChange) > 1,]
interaction_dftop <- interaction_dftop[order(interaction_dftop$log2FoldChange, decreasing=TRUE),]
mat_int <- assay(vsd)[rownames(interaction_dftop), rownames(colData(vsd))]
colnames(mat_int) <- rownames(colData(vsd))
base_mean_int <- rowMeans(mat_int)
mat_int_scaled <- t(apply(mat_int, 1,scale))
colnames(mat_int_scaled) <- colnames(mat_int) # Original code, changed sample names in final figure
num_keep <- 25
rows_keep <- c(seq(1:num_keep), seq((nrow(mat_int_scaled)-num_keep),nrow(mat_int_scaled)))
l2_val <- as.matrix(interaction_dftop[rows_keep,]$log2FoldChange)
colnames(l2_val) <- "logFC"
mean <- as.matrix(interaction_dftop[rows_keep,]$baseMean)
colnames(mean) <- "AveExp"
# maps values between b/w/r for min and max l2 values
col_logFC <- colorRamp2(c(min(l2_val),0, max(l2_val)),c("blue", "white","red"))
# maps between 0% quantile, and 75% quantile of mean values
col_AveExpr <- colorRamp2(c(quantile(mean)[1],quantile(mean)[4]),c("white", "red"))
ha <- HeatmapAnnotation(summary = anno_summary(gp=gpar(fill=2), height =unit(2, "cm")))
h1 <- Heatmap(mat_int_scaled[rows_keep,],cluster_rows = F, column_labels = colnames(mat_int_scaled), name="Z-score", cluster_columns = T)
h2 <- Heatmap(l2_val, row_labels = rownames(interaction_dftop[rows_keep,]), cluster_rows = F, name= "logFC", top_annotation = ha, col= col_logFC, cell_fun = function(j,i,x,y,w,h,col){
grid.text(round(l2_val[i,j],2),x,y)
})
h3 <- Heatmap(mean, row_labels = rownames(interaction_dftop[rows_keep,]), cluster_rows = F, name= "AveExpr", col = col_AveExpr, cell_fun = function(j,i,x,y,w,h,col){
grid.text(round(mean[i,j],2),x,y)
})
h <- h1+h2+h3
h
